
%%% Weighting the Evidence: A Rank-Dependent Model of Outdoor Recreation, June 2024
%%% This file predicts catch outcomes and participation for year 2021 in New Jersey under an increase in the bag limit

%Importing the different data files. 
Tmain=readtable('\ReadyDataFluke2022.xlsx'); 
cost=importdata('trip_costs_NE.xlsx'); 
cost=cost.data(cost.textdata(:,2)=="NJ");
C=unique(cost)/1000; 
y=[Tmain.age Tmain.income_medium Tmain.income_high Tmain.education_college Tmain.education_graduate Tmain.avidity Tmain.male];
Y=unique(y,'rows'); 
CatchPerTrip=importdata('CatchPerTrip2021.mat');
SFreshaffle=CatchPerTrip.SFreshaffle2021;
Sreshaffle=CatchPerTrip.Sreshaffle2021;
BSBreshaffle=CatchPerTrip.BSBreshaffle2021;
clearvars SFreshaffle21 Sreshaffle21 BSBreshaffle21 Tmain
parameters_rdu_tk=importdata('parameters_rdu_tk.txt'); 
parameters_rdu_prelec=importdata('parameters_rdu_prelec.txt'); 
sigma_rdu_tk=importdata('sigma_rdu_tk.txt');  
sigma_rdu_prelec=importdata('sigma_rdu_prelec.txt');
K=importdata('k_NJ_2021.xls');
F_size_SF=importdata('ProbFlukeSizesNJ2021.txt'); 
p_cutoff_SF=importdata('p_cutoff_SF.txt'); 

SFbaglimit=3; 
SFsizelimit=18; 

NETrips=importdata('AvgeHarvestTrip.xlsx');
[r7,c7]=find(contains(NETrips.textdata,'NJ'));
NETrips.data=NETrips.data(r7,:);
[r8,c8]=find(NETrips.data==2021);
NETrips.data=NETrips.data(r8,:);
TRIPS=round(NETrips.data(1,2)/1000);
cutoff_sf=NETrips.data(1,3);
bsb_k=NETrips.data(1,4);
scup_c=NETrips.data(1,5);
sf_r=NETrips.data(1,6);
bsb_r=NETrips.data(1,7);
clearvars NETrips
T=100; 
M=25000; 

[AggregateKeepSFnumber_eu,AggregateReleasedSFnumber_eu,AggregateKeepSnumber_eu,AggregateKeepBSBnumber_eu,AggregateReleasedBSBnumber_eu,AggregateKeepSFnumber_rdu_tk,AggregateKeepSFnumber_rdu_prelec,AggregateReleaseSFnumber_rdu_tk,AggregateReleaseSFnumber_rdu_prelec]=deal(zeros(M,T));  
[AggregateKeepSnumber_rdu_tk,AggregateKeepSnumber_rdu_prelec,AggregateKeepBSBnumber_rdu_tk,AggregateKeepBSBnumber_rdu_prelec,AggregateReleasedBSBnumber_rdu_tk,AggregateReleasedBSBnumber_rdu_prelec,TotalWelfare_eu,TotalWelfare_rdu_tk,TotalWelfare_rdu_prelec,NumberOftrips_eu,NumberOftrips_rdu_tk,NumberOftrips_rdu_prelec]=deal(zeros(M,T));
[vartheta_eu,vartheta_rdu_tk,vartheta_rdu_prelec]=deal(zeros(T,7)); 

tic; 
rng(0, 'combRecursive');
N=30; 

u_eu=zeros(T,N); 
[v_eu,opt_out_rdu_tk,opt_out_rdu_prelec,opt_out_eu,probTripRDU_tk,probTripEU,probTripRDU_prelec,v_rdu_tk,v_rdu_prelec,ReleasedSFnumber,ReleasedSnumber,ReleasedBSBnumber,BSBnumber,Snumber,SFnumber]=deal(zeros(M,T)); 
[SFkept,BSBkept,BSBreleased,SFkeptNumbers,SFreleased,Scatch]=deal(zeros(N,T)); mean_sf=zeros(T,1);

par_rdu_tk=mvnrnd(parameters_rdu_tk,sigma_rdu_tk,T);
par_rdu_prelec=mvnrnd(parameters_rdu_prelec,sigma_rdu_prelec,1.2*T);
par_rdu_prelec(par_rdu_prelec(:,16)<0,:)=[];
par_eu=par_rdu_tk;
par_eu(:,length(parameters_rdu_tk))=1;

nbinR=zeros(T,1);
nbinP=zeros(T,1);
for t=1:T
     x=reshape(SFreshaffle(:,:,t),[],1);
     keep=zeros(1,length(x));
for j=1:length(x)
z=rand(1,x(j));
keep(j)=sum((z>=p_cutoff_SF),'all');
end
pd=fitdist(keep','nbin');
nbinR(t)=pd.R;
nbinP(t)=pd.P;
end

Cost_data=[datasample(C,M) datasample(C,M)];
Indcharac_data=datasample(Y,M,1);

%1. RANK-DEPENDENT UTILITY TVERSKY & KAHNEMAN
  k=0;
  while k<=K.data(1)     
    
    k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end

     end 
    
    vartheta_rdu_tk(t,:)=[par_rdu_tk(t,9),par_rdu_tk(t,10),par_rdu_tk(t,11),par_rdu_tk(t,12),par_rdu_tk(t,13),par_rdu_tk(t,14), par_rdu_tk(t,15)];
    v_rdu_tk(k,t)=RDU_tk(par_rdu_tk(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac);
    opt_out_rdu_tk(k,t)=par_rdu_tk(t,8)+IndvCharac*vartheta_rdu_tk(t,:)'+par_rdu_tk(t,7)*CostCharacOpt_out;
    probTripRDU_tk(k,t)=exp(v_rdu_tk(k,t))/(exp(v_rdu_tk(k,t))+exp(opt_out_rdu_tk(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*SFnumber(k,t);
    AggregateReleaseSFnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*sf_r;
    AggregateKeepBSBnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*bsb_k;
    AggregateReleasedBSBnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*bsb_r;  
    TotalWelfare_rdu_tk(k,t)=log(exp(v_rdu_tk(k,t))+exp(opt_out_rdu_tk(k,t)))/(-par_rdu_tk(t,7)/1000);
    NumberOftrips_rdu_tk(k,t)=probTripRDU_tk(k,t);
    end
  end

varNames1={'Number of trips RDU TK','Number of choice occasions RDU TK','Welfare RDU TK','Number of SF kept RDU TK', 'Number of BSB kept RDU TK'};
writetable(table(sum(NumberOftrips_rdu_tk,1)',k*ones(1,T)',sum(TotalWelfare_rdu_tk,1)',sum(AggregateKeepSFnumber_rdu_tk,1)',sum(AggregateKeepBSBnumber_rdu_tk,1)','VariableNames',varNames1),'PredictedNJ_2021_bl_up.xls','Sheet','RDU_KT_before')
  
%2. RANK-DEPENDENT UTILITY PRELEC
  k=0;
  while k<=K.data(2)   
    
    k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     end 
    
    vartheta_rdu_prelec(t,:)=[par_rdu_prelec(t,9),par_rdu_prelec(t,10),par_rdu_prelec(t,11),par_rdu_prelec(t,12),par_rdu_prelec(t,13),par_rdu_prelec(t,14), par_rdu_prelec(t,15)];
    v_rdu_prelec(k,t)=RDU_prelec(par_rdu_prelec(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac);
    opt_out_rdu_prelec(k,t)=par_rdu_prelec(t,8)+IndvCharac*vartheta_rdu_prelec(t,:)'+par_rdu_prelec(t,7)*CostCharacOpt_out;
    probTripRDU_prelec(k,t)=exp(v_rdu_prelec(k,t))/(exp(v_rdu_prelec(k,t))+exp(opt_out_rdu_prelec(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*SFnumber(k,t);
    AggregateReleaseSFnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*sf_r;
    AggregateKeepBSBnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*bsb_k;
    AggregateReleasedBSBnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*bsb_r;  
    TotalWelfare_rdu_prelec(k,t)=log(exp(v_rdu_prelec(k,t))+exp(opt_out_rdu_prelec(k,t)))/(-par_rdu_prelec(t,7)/1000);
    NumberOftrips_rdu_prelec(k,t)=probTripRDU_prelec(k,t);
    end
  end

varNames2={'Number of trips RDU Prelec','Number of choice occasions RDU Prelec','Welfare  RDU Prelec','Number of SF kept RDU Prelec','Number of BSB kept RDU Prelec'};
writetable(table(sum(NumberOftrips_rdu_prelec,1)',k*ones(1,T)',sum(TotalWelfare_rdu_prelec,1)',sum(AggregateKeepSFnumber_rdu_prelec,1)',sum(AggregateKeepBSBnumber_rdu_prelec,1)','VariableNames',varNames2),'PredictedNJ_2021_bl_up.xls','Sheet','RDU_Prelec_before')

  
%3. EXPECTED UTILITY
    k=0;
  while k<=K.data(3)
        k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     u_eu(t,j)=(1-exp(-par_eu(t,6)*(par_eu(t,2)*SFkept(j,t)/100+par_eu(t,1)*sf_r/100+par_eu(t,4)*bsb_k/100+par_eu(t,3)*bsb_r/100++par_eu(t,5)*scup_c/100)))/par_eu(t,6);   
     end 
   
    vartheta_eu(t,:)=[par_eu(t,9),par_eu(t,10),par_eu(t,11),par_eu(t,12),par_eu(t,13),par_eu(t,14), par_eu(t,15)];

    par_eu=par_rdu_tk;
    z=size(par_rdu_tk,2);
    par_eu(:,z)=1;
    v_eu(k,t)=RDU_tk(par_eu(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac); 
    opt_out_eu(k,t)=par_eu(t,8)+IndvCharac*vartheta_eu(t,:)'+par_eu(t,7)*CostCharacOpt_out;
    probTripEU(k,t)=exp(v_eu(k,t))/(exp(v_eu(k,t))+exp(opt_out_eu(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_eu(k,t)=probTripEU(k,t)*SFnumber(k,t);
    AggregateReleasedSFnumber_eu(k,t)=probTripEU(k,t)*sf_r;
    AggregateKeepBSBnumber_eu(k,t)=probTripEU(k,t)*bsb_k;
    AggregateReleasedBSBnumber_eu(k,t)=probTripEU(k,t)*bsb_r;
    TotalWelfare_eu(k,t)=log(exp(v_eu(k,t))+exp(opt_out_eu(k,t)))/(-par_eu(t,7)/1000);
    NumberOftrips_eu(k,t)=probTripEU(k,t);
    end
  end


varNames3={'Number of trips EU','Number of choice occasions EU','Welfare EU','Number of SF kept EU','Number of BSB kept EU'};
writetable(table(sum(NumberOftrips_eu,1)',k*ones(1,T)',sum(TotalWelfare_eu,1)',sum(AggregateKeepSFnumber_eu,1)',sum(AggregateKeepBSBnumber_eu,1)','VariableNames',varNames3),'PredictedNJ_2021_bl_up.xls','Sheet','EU_before')

SFbaglimit=4; 

%1. RANK-DEPENDENT UTILITY TVERSKY & KAHNEMAN
  k=0;
  while k<=K.data(1)     
    
    k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
    end
    
    vartheta_rdu_tk(t,:)=[par_rdu_tk(t,9),par_rdu_tk(t,10),par_rdu_tk(t,11),par_rdu_tk(t,12),par_rdu_tk(t,13),par_rdu_tk(t,14), par_rdu_tk(t,15)];
    v_rdu_tk(k,t)=RDU_tk(par_rdu_tk(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac);
    opt_out_rdu_tk(k,t)=par_rdu_tk(t,8)+IndvCharac*vartheta_rdu_tk(t,:)'+par_rdu_tk(t,7)*CostCharacOpt_out;
    probTripRDU_tk(k,t)=exp(v_rdu_tk(k,t))/(exp(v_rdu_tk(k,t))+exp(opt_out_rdu_tk(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*SFnumber(k,t);
    AggregateReleaseSFnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*sf_r;
    AggregateKeepBSBnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*bsb_k;
    AggregateReleasedBSBnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*bsb_r;  
    TotalWelfare_rdu_tk(k,t)=log(exp(v_rdu_tk(k,t))+exp(opt_out_rdu_tk(k,t)))/(-par_rdu_tk(t,7)/1000);
    NumberOftrips_rdu_tk(k,t)=probTripRDU_tk(k,t);
    end
  end

varNames4={'Number of trips RDU TK','Number of choice occasions RDU TK','Welfare RDU TK','Number of SF kept RDU TK', 'Number of BSB kept RDU TK'};
writetable(table(sum(NumberOftrips_rdu_tk,1)',k*ones(1,T)',sum(TotalWelfare_rdu_tk,1)',sum(AggregateKeepSFnumber_rdu_tk,1)',sum(AggregateKeepBSBnumber_rdu_tk,1)','VariableNames',varNames4),'PredictedNJ_2021_bl_up.xls','Sheet','RDU_KT_after')
  
%2. RANK-DEPENDENT UTILITY PRELEC
  k=0;
  while k<=K.data(2)   
    
    k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t);
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
     end 
    
    vartheta_rdu_prelec(t,:)=[par_rdu_prelec(t,9),par_rdu_prelec(t,10),par_rdu_prelec(t,11),par_rdu_prelec(t,12),par_rdu_prelec(t,13),par_rdu_prelec(t,14), par_rdu_prelec(t,15)];
    
    v_rdu_prelec(k,t)=RDU_prelec(par_rdu_prelec(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac);
    opt_out_rdu_prelec(k,t)=par_rdu_prelec(t,8)+IndvCharac*vartheta_rdu_prelec(t,:)'+par_rdu_prelec(t,7)*CostCharacOpt_out;
    probTripRDU_prelec(k,t)=exp(v_rdu_prelec(k,t))/(exp(v_rdu_prelec(k,t))+exp(opt_out_rdu_prelec(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*SFnumber(k,t);
    AggregateReleaseSFnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*sf_r;
    AggregateKeepBSBnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*bsb_k;
    AggregateReleasedBSBnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*bsb_r;  
    TotalWelfare_rdu_prelec(k,t)=log(exp(v_rdu_prelec(k,t))+exp(opt_out_rdu_prelec(k,t)))/(-par_rdu_prelec(t,7)/1000);
    NumberOftrips_rdu_prelec(k,t)=probTripRDU_prelec(k,t);
    end
  end

varNames5={'Number of trips RDU Prelec','Number of choice occasions RDU Prelec','Welfare  RDU Prelec','Number of SF kept RDU Prelec','Number of BSB kept RDU Prelec'};
writetable(table(sum(NumberOftrips_rdu_prelec,1)',k*ones(1,T)',sum(TotalWelfare_rdu_prelec,1)',sum(AggregateKeepSFnumber_rdu_prelec,1)',sum(AggregateKeepBSBnumber_rdu_prelec,1)','VariableNames',varNames5),'PredictedNJ_2021_bl_up.xls','Sheet','RDU_Prelec_after')
  

%3. EXPECTED UTILITY
    k=0;
  while k<=K.data(3)
        k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1);
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     u_eu(t,j)=(1-exp(-par_eu(t,6)*(par_eu(t,2)*SFkept(j,t)/100+par_eu(t,1)*sf_r/100+par_eu(t,4)*bsb_k/100+par_eu(t,3)*bsb_r/100++par_eu(t,5)*scup_c/100)))/par_eu(t,6);   
     end 
   
    vartheta_eu(t,:)=[par_eu(t,9),par_eu(t,10),par_eu(t,11),par_eu(t,12),par_eu(t,13),par_eu(t,14), par_eu(t,15)];
    par_eu=par_rdu_tk; 
    z=size(par_rdu_tk,2);
    par_eu(:,z)=1;
    v_eu(k,t)=RDU_tk(par_eu(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac); 
    opt_out_eu(k,t)=par_eu(t,8)+IndvCharac*vartheta_eu(t,:)'+par_eu(t,7)*CostCharacOpt_out;
    probTripEU(k,t)=exp(v_eu(k,t))/(exp(v_eu(k,t))+exp(opt_out_eu(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));

    AggregateKeepSFnumber_eu(k,t)=probTripEU(k,t)*SFnumber(k,t);
    AggregateReleasedSFnumber_eu(k,t)=probTripEU(k,t)*sf_r;
    AggregateKeepBSBnumber_eu(k,t)=probTripEU(k,t)*bsb_k;
    AggregateReleasedBSBnumber_eu(k,t)=probTripEU(k,t)*bsb_r;
    TotalWelfare_eu(k,t)=log(exp(v_eu(k,t))+exp(opt_out_eu(k,t)))/(-par_eu(t,7)/1000); 
    NumberOftrips_eu(k,t)=probTripEU(k,t);
    end
  end

varNames6={'Number of trips EU','Number of choice occasions EU','Welfare EU','Number of SF kept EU','Number of BSB kept EU'};
writetable(table(sum(NumberOftrips_eu,1)',k*ones(1,T)',sum(TotalWelfare_eu,1)',sum(AggregateKeepSFnumber_eu,1)',sum(AggregateKeepBSBnumber_eu,1)','VariableNames',varNames6),'PredictedNJ_2021_bl_up.xls','Sheet','EU_after')
  

toc;
save('PredictedNJ2021_bl_up.mat');



